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ABSTRACT : A subcritical panel method applied to flow analysis and 

aerodynamic design of complex aircraft configurations is presented. 

The analysis method is based on linearized, compressible, subsonic flow 
equations and indirect Dirichlet boundary conditions. Quadratic dipol 
and linear source distribution on flat panels are applied. 

In the case of aerodynamic design the geometry which minimizes differences 
between design and actual pressure distribution is found iteratively using 
numerical optimization technique. Geometry modifications are modelled by 
surface transpiration concept. Constraints in respect to resulting geometry 
can be specified. A number of complex 3-dimensional design examples are 
presented. The software is adopted to personal computers, and as result an 
unexpected low cost of computations is obtained. 

INTRODUCTION 

One of the most important task in aerodynamic design is such airplane shape 
definition which fulfills the following requirements: low CD, high MA dd and 

CL , appropriate boundary layer stability and stall progression, elimination 

M A X 

of shock waves etc. This, however, depends on appropriate pressure 
distribution on the surface. It is extremely difficult to fulfill all these 
requirements for complex, 3-dimensional airplane configurations where strong 
interference effects occur between aerodynamical ly close coupled elements. 
Optimal design of each element does not lead to optimum of configuration 
because of adverse interference effects. But in principle it is possible to 
design such configurations with neutral or even favorable interference, where 
interaction between airplane components gives benefits and leads to better 
global characteristics then those of separated elements. It is impossible to 
realize such a configuration only on the ground of experimental technique. 
Computational methods of aerodynamics, which have developed quickly during 
last 30 years enable, in connection with the aerodynamic concepts worked out 
at this time (“roof-top", "peaky" etc.), to realize many interesting designs. 
The problem can be illustrated by wing-nacelle-pylon configuration. In the 
past the nacelles were shaped as axisymmetrical body and mounted to swept wing 
by plane pylons. A strong adverse interference occurs leading to loss in 
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("isobar sweep, higher local Mach numbers and shock waves, losses in lift 

coefficient at design angle of attack etc., creating the lower aerodynamic 

efficiency. Later the method of designing for neutral interference was worked 

out, where nacelle and pylon were shaped along stream lines of isolated wing 

in order to minimize interference. It is difficult however, even now, to 

design such configurations with favorable interference. 

Slightly simplifying the problem we can consider three kinds of design 

treatments in aerodynamics using computational methods: 

1. Design by trial and error method 

2. Direct optimization method 

3. Inverse design method 

The first is direct transformation of the wind-tunnel technique on the 
computational ground, where wind-tunnel is replaced by computational system 
and the process of "aerodynamic model manufacture" and "testing" is 
significantly cheaper and faster. Experienced aerodynamicist analyses 
results, specifies the needed modifications and the process is repeated until 
satisfactory computational results are obtained. 

In the second method geometry which minimizes aerodynamic object function 
(such as drag) and fulfills additional constraints is found directly without 
external detailed considerations about flow properties. This method, 
conceptually very attractive and fully automated, can not be actually 
performed in the case of complex configurations because of very high cost 
and many times too low accuracy of up-to the date flow analysis methods which 
lead to so called "numerical noise" and make impossible to find real solution. 

The third method is actually the most effective and refined method 
acceptable in practice. It consists of two steps. First is such a pressure 
distribution specification which fulfills aerodynamic requirements. In the 
second step the geometry corresponding to this pressure is calculated using 
inverse method. It is obvious that the possession of the appropriate inverse 
method is worthy. The method presented in the paper is actually probably the 
most general inverse method applied to subsonic flow region, which allows to 
design of real complex configurations even via interference effects. 


FLOW ANALYSIS 

The method is based on linearized theory of compressible flow [1], 
The Prandt 1 -Cl auert equation 
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[is assumed to govern the perturbation velocity potential in the flowfield. 

The linearized mass flux boundary conditions on external surface are applied 

W* n = (V + w)*n = in /p (2) 

00 s 00 

and express the intensity of mass outflow through the surface, 
w is the perturbation mass flux vector defined by 


W = (0 V x .<P y .V z ) 

The second order pressure formula [assuming V = (1,0,0)] 

cp 2 = -2» x - (g 2 » x +»y+y z ) 

is applied to find aerodynamic forces and moments, and isentropic 
formula is used to express pressure distribution on the surface: 

I Z ~r~r » 1 " ~r ^ I 


■{[> * (' - v v ] ' 


Applying Greens Theorem to the flowfield the perturbation 
velocity potential on the surface can be expressed as: 



where <<p> is the jump of potential across the wake and E is function of 
position (respectively: 1, 1/2 and 0 for P in the flowfield, on the surface 

and outside the flowfield). Equation (6) is solved by panel method based on 
quadratic dipol and linear source distribution on flat panels and indirect 
Dirichlet boundary conditions (zero perturbation potential is specified on the 
internal side of surface). Control points and unknown singularity parameters 
are located in panel center of gravity. Jump of potential across the wake is 
determined by Kutta condition: flow behind the trailing edge of lifting 
surface must be tangent to trailing edge bisector. Finally the integral 
equation (6) is replaced by system of linear equations of the form: 
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f which is solved to obtain the perturbation potential on the surface and jump ~j 

of potential across the wake. Velocity distribution on the surface is 

obtained by numerical differentiation of perturbation potential and adding 

the free-stream contribution. In the local panel coordinate system: 


V = d<p/dt + V .t 

t 00 

V = d<p/ds + V .s 

s oo 


INVERSE METHOD 

The inverse problem is solved in the present method via optimization. 

The method is extension of the previous design method of the author. 

The requested geometry of configuration is searched in a form of sum of the 
initial geometry and linear combination of basic design shapes: 


ND 

GEOMETRY = INITIAL GEOMETRY + ^ X. • (i-th BASIC SHAPE) 
i = l 

Coefficients are found from the condition of minimizing 
the error in pressure distribution: 


( 9 ) 


NP 

E ■ I v (cp j- cp 5 |2 

J=i 


(10) 


where: W - weight function of j-th point 
n D . . 

Cp - design pressure coefficient 
j 

Cpj - its actual value 
using numerical optimization technique. 

Direct application of panel method to find the object function brings 
the high cost of computations. In the presented method the basic design 
shapes are modelled by surface transpiration. The mass flux through the 
surface which shift the stream surface with the distance h normal to the 
initial surface is given by: 


1 

' a (pUh) a (pvh) ' 

-f. 


U TR p 

00 

5t) 



The mean value of the transpiration over the panel is obtained by mass flux 
balance in the volume enclosed by body surface and modelled stream surface. 









267 


Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 
(ICIDES — III). Editor: G.S. Puli kravich. Washington D.C. .October 23-25,1991 


r 





The incremental potential distribution due to surface transpiration 

(i-th basic shape) is calculated from linear equations system similar to (7): 


5<f I 



Potential on the new geometry is expressed as the sum of initial potential 
distribution and linear combination of incremental potential distribution due 
to basic shapes. The new velocity distribution is calculated using eq. (7) 
with new potential value and unit tangent vectors taken from the new geometry. 
Geometry redefinition is performed directly using eq. (9). The optimization 
is performed by quadratic programing method. Additionally geometrical 
constraints are introduced via penalty function. Gradient and Hessian of 
object and penalty functions are calculated analytically which lead to high 
accuracy and low cost. Because of nonlinear nature of the design problem it 
is solved iteratively using geometry obtained after actual design iteration as 
initial in the next one. Block diagram of the method is shown on the Fig. 2. 


COMPUTER CODE 

The method described above was coded in FORTRAN 77 language and implemented 
on PC-Computers. Because of hardware limitations it is performed as a 
package of programs. All basic parts of the method are performed by 
separate computer program, which are sequentially started from batch file. 
The software package consists of 13 programs including two methods of 
solution of linear equations system (iterative and block Crout 
decomposition) and post -processing program. The iterative method of 
gojution performs matrix modification and makes possible to use this method 
even when other iterative methods do not provide the convergence. 

It is possible to use up to 1200 body panels, 500 wake panels, 80 Kutta 
points, 1280 unknown singularity parameters (plus symmetry condition), and 
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Fig. 2 Block diagram 


of the design method 


63 basic design shapes. Flow analysis for PC-386/25MHZ and 1000 panels (plus 
symmetry) took about 25 -40 if using iterative method of solution and about 
60’ if Crout decomposition method is used. Design process took about 12’ for 
40 basic design shapes using Crout method. Using computer 486 computing time 
is about 50% shorter. Cost of such computations is unexpectedly low. 


RESULTS 

Flow an alysis . To show efficiency and accuracy of the method results of 
analysis of test cases from AGARD AG-241 are shown on Fig. 4 and 5. Results 
for RAE WING and STRAKED WING with NACA 0002 profile is compared with 
Datum Results of Rubbert and Roberts. 

560 panels were used (40x14) for RAE WING and 640 (40x16) for STRAKED WING. 
Computing time on PC-386/25MHz respectively 10’ (iter)/16.5* (Crout) and 
14’/23’. It is seen excellent agreement with compared methods. 


Full aircraft configuration design . It consists of wing, body, tail and 
rear mounted nacelle and pylon. The geometry of the configuration is shown 
on Fig. 6. A new pressure distribution (of "roof-top" type) is specified on 
the wing upper surface. At all points of pylon where initial negative pressure 
^exceeds Cp = -0.5 this value was specified as design one. 
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|"~ 38 basic design shapes of spline-support type were specified. The idea of 


this type of shapes is shown on Fig. 3. Node lines on the surface in both 
directions are specified and movement of the node of such network in 
specified direction corresponds to the desired shape function. To find 
movement of other points of the surface the interpolation spline is used. 
The shape functions used correspond to: 

-changes of upper surface section of the wing at four control stations 
(wing-body- junction, rj = 0.3, 0.5 and 1.0] corresponding to vertical 
displacement of points with max. laying at 75%, 55%, 40%, 25%, 15%, 9% and 
4% of arc length (measured from leading edge to trailing edge) 

-changes of wing twist at wing-body- junct ion, tj = 0.5 and 1.0 

-changes of fuselage width in the pylon region with max. at four stations 

-changes of nacelle width in the pylon region with max. at three stations 


Fig. 3 The idea 
of spline support 
basic shapes 



Geometrical constraints used: 

-distance between network points near the trailing edge (for control the 
trailing angle) 

-distance between network points near the max thickness (for control the 
thickness ) 

-distance between network points near the leading edge (for control the 
leading edge radius) 

-distance (in vertical direction) between leading edge and trailing edge 
(for control twist) at three control stations 
-distance between points of pylon (at pylon-fuselage intersection) and 
symmetry plane (for control fuselage shape) at three stations 
-distance between points of pylon (at pylon-nacelle intersection) and 
symmetry plane (for control pylon shape) at three stations. 


1042 body panels, 72 wake panels and 1068 unknown singularity parameters for 
half geometry were used. Computing time using PC-386/25: analysis 78’, design 
cycle 14’. Isobar pattern on the initial geometry and after four design 
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[""iterations are shown on Fig. 7. Pressure distribution at four wing sections 

before and after designing are presented on Fig. 8 and pressure distribution on 

the pylon on Fig. 9. The shape of body-nacelle region and isobar pattern is 

shown on Fig. 10. The region of higher negative pressure occurs on the 

fuselage and nacelle in front of pylon. Adding two shape functions modifying 

fuselage in front of the pylon and specifying additional points with design 

pressure, the result (after 4 iterations) as on Fig. 11 can be obtained. 

The convergence history of the design process is shown on Fig. 12. 

Wingz body-underwing nacelle configuration . The geometry of the configuration 
and details of nacelle region are shown on Fig. 13. 1160 body panels, 104 wake 

panels and 1183 unknowns were used. The pressure on the wing-body alone 
configuration was calculated. Results are shown on Fig. 14a (lower and upper 
surface respectively). Pressure distribution obtained for this conf igurat ion 
is used as design pressure for wing-body-nacelle. Adding plane pylon and 
axisymmetrical nacelle the new pressure distribution and isobar pattern are 
obtained: Fig. 14b. Isobar pattern on the wing after four design iterations is 
shown on Fig. 14c, shape of pylon and nacelle on Fig. 15 and pressure at 
subsequent wing sections before and after designing on Fig. 16. Shape of pylon 
section before and after designing is seen on Fig. 17. 38 basic design shapes 

of spline-support type were used. Wing was changed at three control stations: 
T 7 = 0.4, 0.5 and 0.6. Four points on upper surface (x/c=0.03, 0.11, 0.27 and 
0.50) five points on lower surface (x/c = 0.06, 0.17, 0.33, 0.50 and 0.72) and 
twist at each of this stations can vary. Additionally four points of upper 
nacelle contour (x/L = 0.38, 0.50, 0.63 and 0.81) and four points of pylon 
mean line (x/c = 0.25, 0.50 0.75 and 1.00) were changed. The constraints, in 
respect to wing thickness and twist, nacelle shape and pylon modification, 
were specified. It should be noted that despite the constraints used are not 
very restrictive some of them are active. As result, for example, pylon has 
nonzero side force (it had tendency to bend more). Convergence history is 
shown on Fig. 18. It is of value to show some aerodynamic coefficient for the 
configuration: 



Cl , 
wing 

C1 

nac-pyl 

Cl. f i 
total 

^ m total 

wing-body alone 

0. 5098 

- 

0.606 

-0. 1463 

initial 

0.4772 

0.0051 

0.575 

-0. 1482 

designed 

0.5096 

-0. 0008 

0.605 

-0. 1439 


^Computing time using PC-386/25: flow analysis 84*, inverse cycle 15'. 
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p Design of transonic wing . The research wing for jet-trainer type aircraft 
was designed via subcritical equivalent pressure distribution concept [3] by 
the author as a part of research investigations on supercritical wing 
performed at Aviation Institute in Warsaw (unpublished Report of Aviation 
Institute in Warsaw). The supercritical wing section (of slightly 
peaky-type pressure distribution) was designed using finite-difference 
method. Equivalent subcritical pressure distribution for swept wing (sweep 
angle of leading edge 20.7°, at 25% chord 17.3°) was calculated and used as 
design pressure on the upper surface of the wing. The originality of the 
method consist in including the off-design characteristics. By modifying 
constraints it was forced max. pressure peak at high angle of attack and 
low Mach number at about u = 0.4, which suggest separation first at this 
station. If max negative pressure was too high at the station under 
consideration, the higher leading edge radius was enforced by constraints 
(worsening, of course, the pressure distribution) and vice versa. 28 basic 
design shapes were used: five kinds of changes of thickness distribution 
along the chord at five control stations along the span and twist at three 
stations. The geometrical constraints in respect to max thickness, 
trailing edge angle, leading edge radius and twist are utilized. 

480 body panels, 24 wake panels and 492 unknown singularity parameters were 
used. The block diagram of the design process can be introduced as follow: 
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Computing time (386/25): analysis 9’ (Crout), 5* (Iter). In each design 
iteration the flow, at high a, was calculated about 3 times. Resulting isobar 
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I pattern and pressure distribution at three wing sections are presented on 


1 


Fig. 19. Geometrical parameters of the resulting wing are shown on Fig. 20. 
Quite unexpected for swept wing R distribution along the span is seen. Max 
of the leading edge radius occurs at 80% of semispan. Max of pressure peak 
(a=12 , Ma=0.2) occurs at t ) = 0.40. The drag divergence Mach number 
obtained in wind tunnel tests is shown on Fig. 22 and beginning of 
separation on Fig. 23 (unpublished Report of Aviation Institute in Warsaw). 

It is seen good agreement with expectation. 


CONCLUDING REMARKS 

The method presented above shows great versatility in the case of design of 
real, complex configurations. It has nearly no restrictions in respect to 
the complexity of the geometry. The major limitation is the lack of 
possibility to take into account modification of planform of the wing and 
necessity to fix leading edge point (twist can be changed only by moving 
vertically trailing edge point). It is possible to take into account 
interference effects in designing, that allows to obtain specified pressure 
distribution on one element by changing geometry of the other. 

Recently the method has been extended to the case of multi-point optimization: 
the pressure distribution on different parts of the surface can be specified 
for different angles of attack and the design process is performed at once. 

The method is exceptionally cheap and efficient because of implementation 
on PC-computers. The possibility to take into accounts some 
characteristics at off-design conditions via constrains was shown also. 


REFERENCES 

1. Ward G.N. - “Linearized Theory of Steady High-Speed Flow" 

- Cambridge University Press 1955 

2. Kubrynski K. - “A Subsonic Panel Method for Design of 

3-Dimensional Complex Configurations with Specified 
Pressure Distribution" - Proceedings of the III GAMM-Seminar 
in Kiel - Vieweg 1988 

3. Slooff J.W. - "Application of Computational Procedures in 

Aerodynamic Design" - AGARD R-712 1983. 



Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 
flCIDES-llO.Editor:G.S.Dulikravich. Washington D.C. .October 23-25.1991 



Fig. 5 STRAKED-WING/NACA 0002 : comparison of results 
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Fig. 6 
Network geometry 
of the entire 
aircraft 
configuration 


Fig. 7 

Isobar pattern 
on the surface 
of the initial 
configuration 


Fig. 8 

Isobar pattern 
on the surface 
after 4 design 
iterations 
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g. 9 Pressure distribution on the pylon before and after designing 
p = —0.5 specified at all points where this value is exceeded) 


t - tfWO-OODV-tMXLXI-W. 


i rrr -un i r tm 


Fig. 10 Fuselage-pylon— nacelle region before and after designing 


^\\I 


r-wm 






CONVERGENCE HISTORY" 
or WNC— BODY— N>CELLE-TAJL 
OE30N PROCESS 


Fig. 1 1 Fuselage-pylon-nacelle region 
after designing with additional shape 
function modyfing fuselage in front 
of the pylon 


Fig. 12 Convergence 
of the design process 
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Fig. 19 Planform, isobar pattern and pressure distribution on the designed wing 
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Fig. 20 Geometrical parameters 
of the designed wing sections 


Fig. 21 Peak pressure distribution 


IP 1—6 WING i\ 

DRAG DfvtfcSENCC MAC *1 NUMBER 

—0. 1 0 - I M I 1 T I 1 I - r r-TTi n i i 1 i m i i t i i t t tt i > i > i i ' ii tht tt 
0.40 0 50 060 0.10 0.60 0.90 

MA-DD 






M 


’Ml 


:R 


Fig. 22 The drag divergence Mach number Fig. 23 Begining of separati 
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